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We discuss phase transitions and the phase diagram of a classical dimer model with anisotropic in- 
teractions defined on a square lattice. For the attractive region, the perturbation of the orientational 
order parameter introduced by the anisotropy causes the Berezinskii-Kosterlitz-Thouless transitions 
from a dimer-liquid to columnar phases. According to the discussion by Nomura and Okamoto for a 
quantum-spin chain system [J. Phys. A 27, 5773 (1994)], we proffer criteria to determine transition 
points and also universal level-splitting conditions. Subsequently, we perform numerical diagonal- 
ization calculations of the nonsymmetric real transfer matrices up to linear dimension specified by 
L — 20 and determine the global phase diagram. For the repulsive region, we find the boundary 
between the dimer-liquid and the strong repulsion phases. Based on the dispersion relation of the 
one-string motion, which exhibits a two-fold "zero-energy flat band" in the strong repulsion limit, 
we give an intuitive account for the property of the strong repulsion phase. 

PACS numbers: 05.20.-y, 05.50.+q 



I. INTRODUCTION 

In the early 1960s, Kasteleyn [l[ and Temperley and 
Fisher 0, Q studied the classical dimer model (DM) 
defining the statistical-mechanical problem of the cov- 
ering of a lattice by dimers. They treated the DM on, 
for example, the square lattice in the thermodynamic 
limit, and obtained the partition function to give the 
extensive entropy of an ensemble of dimer configura- 
tions. In particular, it was shown that the close-packed 
model defined on planar lattices can be solved exactly 
by Pfaffian techniques 0]. The properties of these en- 
sembles were studied in subsequent research. For ex- 
ample, those defined on bipartite (non-bipartite) lattices 
such as the square (triangular) lattice exhibit critical 
(off-critical) behavior [fl which are now thought to 
reflect an existence (absence) of the height representa- 
tions for the DMs @, @, S, Efl, E3- While the relevance 
of dimers as diatomic molecules adsorbed on a lattice is 
clear and direct, it can be also related to other degrees of 
freedom [l2| • For instance, the zero-temperature Ising- 
spin antifcrromagnet on a triangular (Villain) lattice can 
be related to the DM on a hexagonal (square) lattice, 
where each dimer represents an unsatisfied bond of spins 
0, H, [IH ■ Also widely known is the string representation 
whereby dimer systems under a certain condition can be 
related to loop gases whose configurations are classified 
by winding numbers (see below) 0, H, [13]. This cor- 
respondence has been utilized in discussions of polymer 
systems [l5[. More importantly, Rokhsar and Kivelson 
introduced quantum dimer models (QDMs) to describe 
the valence-bond physics in quantum Hciscnbcrg antifer- 
romagnets, where the dimer represents a tightly binding 
singlet pair of quantum spins [161 ] . 

These are but a fraction of the examples that show 
DMs having importance in a wide range of research and 
drawing attention over the years; in particular, current 
interest has been mainly focused on exotic phases with 
topological orders observed in the QDMs [171 . More re- 



cently, Blunt et al. reported on an adsorption experiment 
of certain rod-like organic molecules on graphite [HI and 
explained its relevance to the DM on a hexagonal lattice 
which includes interactions between neighboring dimers 
[lj|. This exhibits the fact that interaction effects in 
classical dimers arc also important from both theoretical 
and experimental viewpoints (20| . 

In this paper, we investigate an interacting dimer 
model (IDM) defined on a square lattice: suppose the 
lattice constant a = 1 and let A denote the set of lattice 
sites. Then the following reduced Hamiltonian expresses 
interactions between two nearest-neighbor dimers: 

H = - [K h n(k+^,l)n(k + ^,l + l) 

(kd)eA 

+K v n(k,l + hn(k + l,l + -) , (1) 

where lattice sites and lattice bonds are denoted, respec- 
tively, as (k, I) and (k+ i, I) with k, I £ Z. We locate the 
dimer occupation numbers n{k + i,Z) = or 1 (binary 
variables) on the bonds. The first (second) term on the 
right-hand side of Eq. (Q]) represents an interaction be- 
tween parallel horizontal (vertical) dimers [see Fig.QJa)]. 
Defining the local Boltzmann weights as 



exp(K h ), 



(2) 



with subscripts on the if's referring to the corresponding 
dimer pair, the partition function Z(h, v) is then express- 
ible as a summation with respect to the dimer configura- 
tions C on A, 



Z(h,v) 



j2h N ^v N ^ c \ 



(3) 



where Nh(C) [N V (C)] represents the number of plaque- 
ttes with parallel horizontal (vertical) dimers. For large 
h or v, an attractive case, a twofold- or a fourfold- 
degenerate state with columnar order is expected to be 
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I . I A> 1 1 1 1 1 1 
I J I 1 1 1 1 1 1 
II 1 1 1 1 1 1 1 1 1 



FIG. 1: Dimer configurations, (a), (b), (c), and (d) repre- 
sent examples of the liquid, the staggered, the HC, and the 
VC states, respectively. The local Boltzmann weights are 
also given in (a). The string representation for (a) using (b) 
as "reference state" (see text) is given by the gray lines in 
the y direction. In (b), a dotted line in the [11] direction 
indicates a counterclockwise rotation of all dimers along the 
line; a dashed line in the [ll] direction exhibits a clockwise 
rotation. In (c), elements a x and ad of the C4^-point group 
representing reflections about solid lines in the y and the di- 
agonal directions are indicated. 



stabilized [see examples in Figs. QJc) andQJd)]. Mean- 
while, for small h and v, a repulsive case, a highly degen- 
erate phase is stabilized [see an example in Fig. Gib)]. 
For the isotropic case Kh = K v , results of some nu- 
merical calculations are already available in the litera- 
ture [13, [ll|, but its extension to the anisotropic case 
Kh K v is still lacking. Therefore, we shall clarify the 
global phase diagram and provide evidence corroborat- 
ing the properties of the phase transitions observed in 
anisotropically interacting dimers. 

According to the effective field theory discussed by 
Papanikolaou, Luijten, and Fradkin (PLF) [2^, one 
effect absent in the isotropic system is a perturba- 
tion by the orientational order parameter; this brings 
about Berezinskii-Kosterlitz-Thouless (BKT) transitions 
[H, HIl- Another effect is a renormalization of the 
so-called geometric factor, which becomes important 
in numerical calculations of universal quantities, e.g., 
the central charge and the scaling dimensions of oper- 
ators (see below). In such cases, as demonstrated in 
our own research on an antiferromagnetic Potts model 
with anisotropic next-nearest-neighbor couplings, the so- 
called level-spectroscopy analysis [ID, [26| can provide an 
effective way to determine phase transition points [27| . 
For this reason, we shall also employ the same strategy 
for the present model [28| . 



For later convenience, we shall briefly explain here the 
string representation of the DM on A. As explicitly ex- 
plained in Ref . [29| , the transformation of a dimer config- 
uration, e.g., Fig. [TJa), to a string configuration is per- 
formed via an XOR operation with reference configura- 
tion [Fig. QJb)]. The XOR operation takes the exclusive 
OR between occupation (binary) numbers in these two 
configurations over each bond. Consequently, we obtain 
strings running in the y direction [see the two gray lines in 
Fig. QJa)]. Due to the close-packing condition, they have 
no end points, and thus the string configurations for a 
Lx L system (L is an even number) with periodic bound- 
ary conditions can be characterized by winding numbers 
(N x ,N y ) satisfying < N x , y < L. While in the numer- 
ical calculation of the transfer matrices we shall employ 
N y as a conserved quantity in the row-to-row transfer of 
configurations, we would rather use a quantity 



M = N y — L/2 



(4) 



(\M\ < L/2) for convenience in our discussion. 

The organization of this paper is as follows: In Sec. HII 
we review previous research results to give an effective 
description of the low-energy and long-distance behavior 
of the IDMs 22j. In particular, the operator content of 
the theory, including expressions of local order parame- 
ters and defect operators, and their scaling dimensions 
are explained in detail. We then provide the conditions 
to determine the BKT-transition points and clarify some 
universal relations among excitation levels, which serve 
as a check of our calculational results. In doing this, 
a correspondence with a frustrated quantum-spin chain 
system plays a guiding role. Thus, this correspondence 
will be emphasized and referred to when appropriate. In 
Sec. IIII1 we summarize our numerical study and results 
obtained by the transfer-matrix calculations based on the 
conformal field theory (CFT) [3(| ■ First, we demonstrate 
that the theoretical predictions in Sec.|TT]can be observed 
precisely via numerical analysis of the excitation spectra. 
Next, we provide the global phase diagram of interacting 
dimers, which includes the BKT, the second-order, and 
the first-order transition lines. Also, in a strong repulsion 
region, we expect a highly degenerate phase including the 
staggered state. We calculate the string-number depen- 
dence of the free-energy density for the isotropic case. 
Furthermore, we investigate the "dispersion relation" of 
the one-string motion, and then based on these data we 
shall try to give an insight into properties of the strong 
repulsion phase. The last section, Sec. IIV[ is devoted to 
discussion and summary. We also provide our method 
to evaluate the original BKT transition in the isotropic 
system. Finally, we compare our data with previous re- 
search results. 



II. THEORY 

Continuum field theories offer unified approaches to in- 
vestigate phase transitions in interacting systems on lat- 
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tices. They are derived in the scaling limit, a — > while 
keeping x = (x±, x 2 ) = (ak, al/Q finite. Here £ is the ge- 
ometric factor, taking a fixed value, e.g., £ = 1 (2/V3) for 
isotropic systems on a square (triangular) lattice. How- 
ever, for anisotropic systems, renormalization of Q is nec- 
essary due to interactions leading to non-universal val- 
ues. The renormalized C, can be also related to the veloc- 
ity of an elementary excitation observed in Tomonaga- 
Luttinger liquids |3l||. Thus, ( disappears from the theo- 
retical description if we properly employ its renormalized 
value; but as we will see in Sec. it becomes rather im- 
portant in numerical calculations. 

According to PLF, the effective description of the IDM 
takes the form of a sine-Gordon field theory [2(j. In the 
present case, its expression is given by the Lagrangian 
density £ = £q + £2 + £4 with 



r _ V2 
1-2 



£i 



: cos2V20 :, 
: cos4%/20 : 



(5) 
(6) 
(7) 



We denote the course-grained height field in the two- 
dimensional Euclidean space as </>(x), which satisfies a 
periodicity in height space of V20 = \Z2(f> + 2nN (N G Z) 
Q. Due to the close-packing of the dimers, the defect 
operators given in terms of the disorder field 9 dual to 
</> are absent from £, so that it represents a roughening 
phase or flat phases of an interface model in three di- 
mensions. The theoretical parameter K (the Gaussian 
coupling) describes the stiffness, and determines the di- 
mensions of the operators on the Gaussian fixed line £o- 
In our notation, the vertex operator with m electric and 



n magnetic charges is given by O n 
whose dimension is 



Kn 



(8) 



Therefore, K = 1 (K = 4) represents the condition that 
the perturbation £2 {£4) becomes marginal. 

To make our discussion more concrete, we introduce 
the average (a) and the difference (d) of the couplings as 

K a , d = I (K v ± K h ) 



where the first (second) subscript refers to the upper 
(lower) sign. Then, around the non-interacting point, 
the parameters in £ are roughly given by 



K 



c x K a: y 2 ~c 2 K d , and y 4 



-C3 



(10) 



(ci,2,3 > 0)- We see that the attractive interaction 
K a > increases K, and tends to stabilize the columnar 
states. Also, 1)2 in £ 2 (the oricntational order parameter) 
is proportional to the difference, while y 4 in £4 which 
is a remnant from the discreteness of the square lattice 
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FIG. 2: (Color online) A schematic representation of the 
BKT RG-flow diagram around the origin of the (yo, 2/2) plane. 
The coordinate frame of the average of couplings K a and the 
difference of couplings (see text) is present as an inset. 



is almost constant. For K > 4, both nonlinear terms are 
relevant, but they are not competing against each other, 
so the fourfold-degenerate columnar state stabilized by 
£4 is only lifted to realize twofold-degenerate columnar 
states by £ 2 (see below) . Since the BKT transition by £4 
was already discussed in the literature [13, Hl[ , we shall 
focus our attention on the role of £2- 

According to the standard argument [12], the 
renormalization-group (RG) flow diagram of the sine- 
Gordon model £0 + ^2 (here £4 is irrelevant) is expressed 
by the BKT RG equations (23|, [24[ ; we depict it by em- 
ploying the coupling constants yo (= 2 — 2K ) and ?/2 
in Fig. [2l where the separatrices y 2 = TUo separate the 
dimer-liquid phase from two types of twofold-degenerate 
columnar phases, namely, the horizontal columnar (HC) 
state consisting of dimers in the horizontal direction and 
the vertical columnar (VC) state consisting of dimers in 
the vertical direction. Now, we can point out that our 
task to treat the BKT transitions in the IDM can be re- 
(9) lated to the investigation of the spin-i XXZ chain with 
next-nearest-neighbor interaction because these share the 
same effective description [2^, [2(|. To specify the rela- 
tionship, we introduce the following operators: 



O = V2 cos V2cj), 
0i,2 = cxp(±i v / 26»), 
3 = v^sinV^. 



(11) 
(12) 

(13) 



Here Oo,3 stand for the horizontal and the vertical com- 
ponents of the columnar local order parameter, and take 
expectation values (Oo) 7^ and (O3) = ((Oo) = 
and (O3) ^ 0) in the HC (VC) phase; Oi, 2 are the de- 
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TABLE I: Discrete symmetries of the principal operators 
(O4.5 are referred to in Sec. IIV|I . The expressions M, k x , 
and P represent the string number, the momentum in the 
x direction, and the parity for the reflection a x , respectively 
(see text). Also given are notations and identifications in both 
dimer and quantum-spin (optional) languages. 



Notations 


Operators 


Identifications 


M 




P 


O 


^cos^<^ 


HC (Neel) 





7T 


-1 


Ol,2 


exp(±iV26) 


monomer (doublet) 


±1 





+ 1 


O3 


V%sin-\/20 


VC (dimer) 








+ 1 


O4 


V2 cos 2y/2(j) 


Orientational 








+ 1 


C>5 


s/2 sin 2y/2(j) 


Plaquette 





7T 


-1 



feet (or monomer) operators which change the winding 
numbers classifying dimer configurations. Alternatively, 
in the quantum-spin chain language, Eqs. (jlip - (|13p cor- 
respond to the Neel, the doublet, and the dimer oper- 
ators, and give the lowest excitations in the Neel, the 
spin-liquid, and the dimer phases, respectively (see Ta- 
ble Hj). Nomura and Okamoto (NO) provided criteria to 
determine the BKT-transition points between the spin- 
liquid and the Neel (or dimer) phases based on one-loop 
calculations of the scaling dimensions of these operators 
[33l ]. Therefore, following their argument, we shall dis- 
cuss procedures to determine the BKT-transition points 
in our IDM. 

Consider a system with a finite- strip geometry, viz., a 
narrow band of width L along the x direction and infi- 
nite length along the y direction. The periodic bound- 
ary condition is imposed across the width of the strip. 
The finite-size corrections to the scaling dimensions of 
the above operators are our key quantities to be evalu- 
ated analytically and numerically. Here, we first consider 
the system near the separatrix 7/2 = — yo where a small 
parameter t can be introduced, so that 7/2 = — yo(l + t) 
(\t\ <C 1). Next, the conformal perturbation calculations 
of the renormalized scaling dimensions were performed 
using the sine-Gordon Lagrangian density, for which the 
results can be summarized as follows: 

x °-l~ \vo{l) (1 + 2t ) > ( 14 ) 

xi,2 - 5 ~ \yo( l )> ( 15 ) 

*3 - \ + \yo(l) (3 + 2t) , (16) 

(I = In L is the logarithmic scale length) [33| . According 
to the discussion by NO, we can find the criterion to 
determine the BKT-transition point t = (i.e., the level- 
crossing condition) 

Xq = Xi )2 (17) 

and the level-splitting condition as 

3zo,i,2 + x 3 _ 1 

4 ~ 2' 1 ' 



The latter is one of the universal relations among the 
excitation levels on the separatrix [34l |. and enables us 
to check the consistency of the calculations. Second, we 
investigate the system near the separatrix yi = yo; it 
proceeds in an analogous way to the above. Writing yi = 
yo(l + 1)> w e then obtain the dimensions as 

x ~ l + l yo Q)(3 + 2t), (19) 

H, 2 ~ l - - ±yo(l), (20) 

x 3 ~ 1-1^(0(1+2*). (21) 

Thus, the level-crossing condition needed to determine 
the transition point is given by 

Xi,2 =x 3 , (22) 

and the level-splitting condition is given by 

Xo+3 ^ 3 = 1. (23) 
4 2 17 

In analogy to the quantum-spin chain, each level crossing 
represents an emergence of a SU(2) multiplct structure 
consisting of the singlet and the triplet states (e.g., X3 
and £0,1,2 at 2/2 = —yo)- Since these are the low-energy 
levels in the level-1 SU(2) Wess-Zumino-Witten model 
[13], our criteria, Eqs. (jTTJ) and (f2"2")l . are natural and also 
convincing from this viewpoint. 

At this stage, two comments are in order about the ad- 
vantage in using these relations (the level-spectroscopy 
approach) and the structure of the phase diagram. In 
Sec. Ill, we will outline the numerical transfer-matrix 
calculations performed to obtain the phase diagram. Al- 
though it can treat systems with a strip geometry, ac- 
cessible sizes are strongly restricted to small values, e.g., 
L < 20 in our calculations. In the BKT transition, as 
seen above, the correction terms are typically given by 
the logarithmic form yo ~ 1/ ln(L/Lo)- If we employ the 
standard KT criterion such as xq = \ to determine the 
transition point, its finite-size estimates include these, 
and thus exhibit a slow convergence in their extrapolation 
to the thermodynamical limit. Alternatively, criteria p7|) 
and ([22]) take the logarithmic corrections into account, so 
they provide finite-size estimates with fast convergences 
[1^1 . Consequently, we can employ the following least- 
squares-fitting form in extrapolating the finite-size data 
Q(L) to the thermodynamic limit L — > 00: 

Q(L) ~ Q(oo) + a/L 2 + b/L\ (24) 

which includes the 1/L 2 term stemming from the x = 4 
irrelevant operators as the leading universal correction 

The biases from the finite-strip geometry disappear in 
the limit, and the phase diagram is symmetric with re- 
spect to the isotropic line Kh = K v which is one of the 
inherent properties of the model. Here, we describe how 
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the symmetry of the lattice model is embedded in the 
sine-Gordon field theory. We shall consider the genera- 
tors of the C4„-point group: the tt/2 rotation (C4) and 
the reflection in the x axis (a x ) about the original site 
[see Fig. QJc)]. As well as coordinate transformations, 
these bring about the following changes in the height field 

C 4 : V20 -> y/2cj> - tt/2, a x : V2(f> -> tt - y/2<j>. (25) 

All other elements are obtained from Eq. (|25[k among 
them, we shall investigate transformations of the La- 
grangian density and the principal operators by reflection 
about the diagonal line, ad (= a x o C4) [see Fig. [He)], 
which shifts the field as 



V2(f) -> tt/2 - V2(j). 



(26) 



Since the orientational order parameter Li is odd for ad, 
the Lagrangian density transforms as 



a d : C(K,y 2 ,yi) ->• C(K,-y 2 ,yi) 



(27) 



This indicates a connection between the positive and the 
negative values of Kd- In addition, the above four oper- 
ators transform as 



O1.2, o 3 ^o . 



(28) 



Thus, the symmetry operation ad interchanges the roles 
of the HC and the VC operators while leaving unchanged 
the doublet of the monomer excitations. Consequently, 
as expected, the level-crossing and the level-splitting con- 
ditions for Kd < 0, Eqs. (|17|) and l|18p. arc translated to 
those for K d > 0, Eqs. J22J and J23J). In Sec. Hill we shall 
provide some numerical data to check this symmetry. 



III. NUMERICAL CALCULATIONS 

Now, consider a system on A with the L x 00 stripe 
geometry and introduce the transfer matrix Tjy/ (L) con- 
nccting nearest-neighbor rows in the y direction [see Fig. 
(Ha)] [13, As mentioned in Scc.Hl the string number 
in the y direction, M, is a conserved quantity, which can 
thus be specified explicitly. We denote the eigenvalues 
as X P (L) and their logarithms as E P (L) = —ln\X p (L)\ 
(p specifics an excitation level such as those listed in the 
above). Then, the conformal invariance provides direct 
expressions of the central charge c and the scaling dimen- 
sion x p in the critical systems as (36l . |37| 

7T 2lT 

E g (L)c±Lf-—c, AE p (L)~—x p . (29) 

Here, E S (L), AE p {L) [= E p (L) - E e (L)], and / corre- 
spond to the ground-state energy, an excitation gap, and 
a free-energy density, respectively. The ground state is 
found in the M = (N y = L/2) sector, and the ex- 
cited levels are also in the sectors specified by the dis- 
crete symmetries given in Table [U In addition, since, 



independent of the value of K, the scaling dimension 
of a level- 1 descendant is equal to 1, it has been uti- 
lized to estimate a velocity of elementary excitation in 
the Tomonaga-Luttinger liquid (see, for example, [25[). 
According to their treatment, the effective geometric fac- 
tor (i.e., inverse velocity) can be also calculated from the 
descendant level, say Eq, as (2?| 



C 1 = lim 



2tt/L 



(30) 



The corresponding excitation with small momentum can 
be found numerically. In calculating c and x from the 
excitation gaps via Eq. (|29p . an estimate of ( first needs 
to be obtained. However, it is not necessary in deter- 
mining the BKT-transition points by Eqs. (JTTJ) and (|22[) . 
because these are homogeneous equations of x, and thus 
the gaps — instead of dimensions — can be used. This is 
one of the advantages of the level-spectroscopy approach 

In the following, we shall provide our results from 
numerical calculations for systems up to size L = 20. 
The methodological aspects of transfer-matrix calcula- 
tions have been well explained in the literature [2l| . Fur- 
thermore, due to the sparse nature of the matrices, we 
can output all elements to hard disk. Then, using the 
ARPACK library [38[ , we can calculate the dominant eigen- 
values of the nonsymmetric real matrices. 

As a demonstration we give the Kd dependence of the 
scaling dimensions at K a = in Fig.[3]Ja), where xq, x\p, 
and X3 are plotted by solid, dotted, and dashed lines, re- 
spectively. While the data are for a system with L = 20, 
we can see the excitation spectra approach quite close to 
the exact ones at the non-interacting point Kd = 00]. 
Furthermore, we find that the HC and the VC excitations 
interchange their behaviors at Kd = 0, and the former 
(the latter) shows a level crossing with the doublet exci- 
tations at a certain negative (positive) value. According 
to theoretical predictions (JTTJ) and ([2"2"| . these can provide 
finite-size estimates of the BKT-transition points to the 
HC and the VC phases; we shall give some evidence to 
support our augment. In Fig. EJb), we exhibit extrapo- 
lations of the finite-size estimates to the thermodynamic 
limit according to Eq. (|24p . The downward (upward) 
pointing triangles exhibit values of —Kd (Kd) at which 
the crossings between xq and x\ y <z (^1,2 and X3) occur in 
the systems with L = 16, 18, and 20. Their extrapolated 
values strongly agree with each other (i.e., their devia- 
tion is within 0.01%), which is the obvious condition to 
be satisfied. In Figs. [3Jc) and E^d), we plot averaged 
values, i.e., the left-hand sides (LHSs) of Eqs. (fT5)l and 
([23| . as well as the dimensions at the BKT-transition 
points estimated in Fig.[3Jb). In both cases, the extrap- 
olated values of the averages agree with the theoretical 
value of |, which exhibits universal level splittings due to 
logarithmic corrections expected at the transition points. 
These observations show that both our strategy and nu- 
merical procedure are valid also for investigations of the 

idm da. 
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FIG. 3: (a) An example of the Kd dependences of scaling 
dimensions at K a = (for L = 20 case). An inserted key 
identifies excitations and lines, (b) The finite-size estimates 
of the transition points to the HC (VC) phase denoted by 
downward (upward) pointing triangles are extrapolated ac- 
cording to Eq. (|24|l (see solid lines), (c) and (d) Checks of 
the universal level-splitting conditions (Tig]) and ([23]) at the 
BKT-transition points. An inserted key identifies excitations 
and marks; x av in (c) [(d)] is the LHS of Eq. (JTSj [Eq. lf23])]. 
and the least-squares fit solid line is exhibited. 



In Fig. [4l we summarize our results of numerical cal- 
culations for the global phase diagram of our model (JTJ) , 
in which the diagonal line (the center point) corresponds 
to the isotropic (noninteracting) system. Due to the fact 
that the phase diagram is symmetric about the line, it 
is sufficient to explicitly calculate one side of the whole 
parameter space; our calculations are thus restricted to 
the region Kd < (i.e., the upper-left triangular area). 
The open circles with solid lines give the phase bound- 
aries between the dimer- liquid and the columnar phases. 
In the area apart from the BKT-transition boundaries, 
we can estimate the Gaussian coupling from the rela- 
tion K = X1.2/X3. The contour lines of A' = ^, |, 
and in addition to the points (plus marks) associated 
with K = 1, 2, and 3 on the isotropic line are given in 
the figure. As expected, the dimer-liquid phase spreads 
over the area satisfying the condition < K < 4. For 
1 < K < 4, this phase only survives on the isotropic line, 
but it is eventually terminated by another BKT transi- 
tion caused by the marginally relevant £4 perturbation. 
We provide our estimation of the point by our approach 
(double circle in the figure) althou gh s ome numerical re- 
sults were previously available [20I l2l| . We shall explain 
our method and compare our results with these in Sec. 
llYl 

To check the criticality of the dimer-liquid phase, we 
estimate the central charge along the line Kh — by the 
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FIG. 4: (Color online) Global phase diagram. The diago- 
nal line (the center point) corresponds to the isotropic (non- 
interacting) system. The circles with solid lines separate the 
dimer-liquid phase from the HC and the VC phases. The 
squares with solid lines exhibit the condition K = 0, which 
is the boundary of the dimer-liquid phase. Also plotted are 
the contour lines of K = g, |, and | as well as the points 
(plus marks) of K = 1, 2, and 3. The double circle indicates 
another BKT-transition point brought about by £4. 

use of relations (|29|) and (f30|) . In Fig. El we give the 
K v dependencies of the effective geometric factor ( (di- 
amonds), the coefficient of the 1/L correction 7 (= c/Q) 
(squares), and their product to estimate c (circles). With 
increasing anisotropy, £ deviates from the isotropic value 
of 1 and approaches a certain value around 1.5 in the limit 
K v — ► — 00. Simultaneously, 7 declines in value and thus 
cannot itself give the universal amplitude of the finite-size 
correction. However, as expected, their product main- 
tains a value c = 1 within the dimer-liquid phase, and 
hence the proper normalization using the effective geo- 
metric factor is necessary for anisotropic systems. In the 
attractive region, one finds a point at which the central 
charge exhibits a steep decrease. We can check that the 
point is almost on the phase transition boundary to the 
VC phase (see the vertical arrow), and thus that it is 
consistent with the level-crossing calculations. 

The dimer-liquid region with 1 < K < 4 corresponds 
to the unstable Gaussian fixed line; the £2 perturbation, 
except for K = 1 , brings about second-order phase tran- 
sitions to the columnar phases. In this case, as l/£ oc 
lA'dl^, the critical exponent characterizing the diverging 
correlation length is given by 1/v = 2 — X2.0 = 2 — 2/K 
[39l ] . To treat this transition, we have performed a finite- 
size-scaling analysis of the corresponding excitation gaps 
and have checked that the scaling behavior is very good 
although we do not provide the data here. In contrast, 



7 



0.5 











■ 1 










o 








<> 










o 


o : c 


























o 






























o 










J w \J U \J »■* 






s 




D □ 










□ 

□ 

□ 






□ 




□ 










□ 






o 






□ 

□ 




o 






□ 












1 












(1-v)/(1+v) 

FIG. 5: (Color online) The K v dependence of £, 7, and c at 
Kh = 0. An inset identifying marks and physical quantities 
is given. The phase boundary between the dimer-liquid and 
the VC phases is around K v ~ 0.62f (see the arrow). 




FIG. 6: (Color online) The free-energy density /m as a func- 
tion of the string number M and the interaction K a for the 
isotropic system with L = 20. The phase boundary between 
the dimer-liquid and the strong repulsion phases is estimated 
as K a ~ -1.97. 



the liquid phase is absent in the more attractive region, 
whence the phase transition between the HC and the VC 
phases becomes first order (the solid line) accompanied 
by a jump in the phase-locking point (y2<f>) from or 7r 
to tt/2 or 3tt/2. 

Finally, we discuss the transition to the strong repul- 
sion phase (the upper-right gray-color region) at which 
the stiffness of the Gaussian model vanishes (see squares 
with solid lines in Fig. 2]). In terms of the height model, 
this vanishin g pe rmits the interface to tilt globally with- 
out cost pfil. |40|| . To get some deeper insight, we shall 
here focus again on the analogy to a transition observed 
in the quantum-spin chain. The spin-i XXZ chain is 
solvable [U, [42| and exhibits c = 1 criticality for the 
anisotropy parameter satisfying — 1 < A < 1 . This phase 
is terminated at the SU(2) ferromagnetic point A = 1, 
where there occurs a first-order phase transition accom- 
panied by the vanishing of the Gaussian coupling. At 
this point, the ground state of the L-site chain forms 
a SU(2) multiplct with total spin L/2 and thus pos- 
sesses L + 1 degeneracy with respect to the z compo- 
nent S% otsl G [—L/2, L/2]. This degeneracy is also im- 
plied from the following theoretical observation: since 
the vertex operator corresponding to 0o,n in Sec. |TT] ex- 
presses a n-spin flip excitation from the ground state with 
S? ota i = 0, the charge n in a L-site system is restricted 
to values in [—L/2, L/2]. In addition, since the scaling 
dimension Xq^ becomes zero in the limit of K — > 0, at 
least the corresponding L excited levels should degener- 
ate to the ground state to realize the L + 1 degeneracy 
[43l . I4II ]. Returning to our DM where the string num- 
ber plays a role as a total magnetization, the ground- 
state energy in each topological sector, Em, s (L), is ex- 
pected to become independent of M. To see this degen- 
eracy, we calculate the M -dependent free-energy density 
L/m(L) = Em, s {L). In Fig. [SI we give the results for the 



isotropic system with L = 20. The average of couplings 
varies within — 3 < K a < and our estimate of the tran- 
sition point is K a ~ —1.97 (see Sec. II Vp . One can see 
that, with decreasing K a , the free energy tends to show a 
weaker M dependence and, for K a smaller than roughly 
the transition point, it becomes almost constant and zero. 
The above-mentioned O(L) degeneracy inferred from the 
instability of the Gaussian criticality seems to be con- 
sistent with this M independence. However, there still 
exists a discrepancy in the degree of degeneracy with the 
staggered state; we shall discuss this issue for the rest of 
this section. 

It is known that the degeneracy of the staggered state 
is subextensive, i.e., tx exp(aL) [2l|, |4f|. This is be- 
cause, as depicted in Fig. [IJb), the tt/2 counterclockwise 
simultaneous rotation of all dimers along a dotted line 
in the [11] direction can be performed independently of 
each other — the same holds also for clockwise rotations 
along the lines in the [11] direction (a dashed line is an 
example). Thus, if one chooses a certain direction, the 
staggered state is completely ordered in that direction 
and completely disordered in the other one. While the 
problem of how the staggered state is stabilized is quite 
unclear, we shall try to give an insight based on an analy- 
sis of the one-string motion. For convenience, we consider 
the transfer matrix connecting the next-nearest-neighbor 
rows, i.e., Tf_ L ^ 2 , so its eigenvalues or their logarithms 

are squared (i.e., A 2 ) or doubled (i.e., 2E), respectively. 
We treat two sites as one unit in which four states are in- 
cluded. Then, we can analytically diagonalize the matrix 
by a Fourier transformation and obtain the g-dependence 
of the energy 2E q , i.e., the "dispersion relation" of the 
one-string motion in the x direction. Here, we only show 
results; the details of how to construct the transfer matrix 
and also the calculation of eigenvalues in the one-string 
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FIG. 7: Analytical results for the dispersion relations of a 
one-string motion. The lower two of the four bands are drawn 
for the isotropic systems. The complex values are two-fold 
degenerate and denoted by dotted lines; the flat band appears 
in the strong repulsion limit [panel (d)]. 



sector are given in the Appendix. In Fig. for sev- 
eral values of the interactions (isotropic cases), we draw 
the lower two of four bands. When the eigenvalue be- 
comes a complex number, we take its magnitude in the 
plot. While the symmetric two-band structure for the 
non-interacting case [Fig.[JJa)] [H,[2(| is deformed by in- 
teractions, there is a unique minimum at q — for finite 
interaction cases [Figs.[7Jb) and[JJc)]. At the same time, 
as expressed by the dotted lines, complex-conjugate pairs 
of eigenvalues start to appear near the zone boundary 
points q = ±7r. And then, in the strong repulsion limit 
[Fig-Hid)], we find an emergence of a two- fold degenerate 
zero-energy flat band. The corresponding eigenvalues are 
given by \ 2 q = e ±lq , which represents the modulations 
with wave number ±q in the y direction [46| . Conse- 
quently, our results show dispcrsionlcss motion in the [11] 
and the [11] directions, which precisely reflect the above- 
mentioned degeneracy, and thus this flat band structure 
may be a signature of the subextensive degeneracy in the 
staggered state. If we accept this naive argument, we can 
conjecture that the staggered state is only realized in the 
limit. However, our argument is of course at a very spec- 
ulative level; a full understanding should include also the 
degeneracy in many-string sectors. 



IV. DISCUSSIONS AND SUMMARY 

The isotropic case was discussed in detail in Refs. [2(| 
l2lT ] , where the BKT-transition point driven by £4 and the 
first-order transition point to the strong repulsion phase 
were numerically obtained. We have also estimated these 



transition points (see the double circle and the double 
square in Fig. 5]); in particular, for the former, the above 
level-spectroscopy approach has been applied. Thus, here 
we briefly explain our procedure and compare the results. 
Since the Lagrangian density Co + £4 is analyzed in the 
region K ~ 4, we focus our attention not on 0o,3 but 
instead on the following order parameters responsible for 
the breaking of the tt/2 rotational symmetry: 



4 = V2 cos 

5 = \/2sin2v / 2(/). 



(31) 
(32) 



While the former is the oricntational order parameter, 
the latter represents the plaquette order 47| which has 
not been found in classical DMs 0, [21 1 (the locking 
points are (y/2(j>) = n/4, 37r/4, 57r/4, and 7n/4). One 
then finds that via the transformation 2\/2<j) — > \[2<\> an d 
K/4 — ► K , the Lagrangian density and the operators are 
reduced to 



C + £4(2/4) -> £0 + £2(2/4), 4 , 5 -> O ,3- 



(33) 



Therefore, from the discussion in Sec. [Ill the level- 
crossing and the level-splitting conditions (fTT)) and l|T8j) 
are satisfied by the scaling dimensions of these operators, 
say Xi t 5 (here, we have taken the condition j/4 < into 
account). Since the half-charge excitations exp(±ji \/29) 
are absent in our system, we employ the condition 



3.T4 + :r,5 



(34) 



to determine the BKT-transition point. The correspond- 
ing excitation levels can be found in the sectors specified 
by their symmetries given in Table |TJ We extrapolate 
the finite-size estimates up to L = 20 to the thermo- 
dynamic limit according to Eq. (|24[) . We then obtain 
the BKT-transition point as K a ~ 1.523. As we see 
in Table HU the agreement with previous results is very 
good, which indicates that our approach is valid. Simi- 
larly, for first-order phase transition, we have estimated 
the point via condition K = (see also Ref. [2l[) while 
others have determined this transition from a point of 
breakdown in the condition c = 1. Our result is closer 
to the estimate of Alet et al. [2(| although there still 
exists considerable discrepancy among these estimates. 
Likewise for instance for phase-separation transitions ob- 
served in one-dimensional electron systems, higher-order 
corrections have been argued to ambiguously affect es- 
timations [44j|. Thus, we think that the discrepancy in 
these estimates may reflect their effects. 

To summarize, we investigated the anisotropically in- 
teracting dimer model on a square lattice. For the at- 
tractive case, the orientational-order-parameter pertur- 
bation introduced by the anisotropy brings about the 
BKT transition to the columnar phases. We pointed 
out the close relationship of our model to a frustrated 
quantum-spin chain and then found the criteria to deter- 
mine the transition points. Using these, we performed 
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FIG. 8: (Color online) The string vacuum |0) and the one-string states {|^4z}, \Bi), \Ci), \Di)} are depicted in the first line. A 
square given by a dotted blue line indicates a unit cell which includes two sites and four bonds. Dimers and strings are given 
by black rectangles and gray lines, respectively. In subsequent lines, 24 microscopic processes of transfers of one-string states 
between two next-nearest-neighbor rows are given with weights. 



TABLE II: Estimations of the BKT and the first-order tran- 
sition points in the isotropic system. For the BKT transition 
the rotational order parameters were treated in Ref. [2(J. In 
Refs. [2jj, HlJ , the first-order transition point were estimated 
from the breakdown of the condition c = 1. 



Criteria 



BKT Criteria First-order 



Ref. [20] Order parameters 1.54 c / 1 -2.23 
Ref. [21] c, x, etc. 1.5-1.7 c/1 -1.39 

Present Equation (J34} 1-523 K = -1.97 



phase boundary, there exist some points with unclear 
status within the strong repulsion phase including the 
staggered state. Based on the dispersion relation of the 
one-string motion, we gave a possible scenario for the 
stabilization of the staggered phase. However, although 
this issue still remains an open question, we now think 
that the nature of the nonsymmetric real matrix might 
have relevance to its description |46| . 



level-spectroscopy analysis of the eigenvalue structures 
of the transfer matrices. Numerical results were then 
summarized as the global phase diagram (Fig. [4j), which 
includes the dimer-liquid, the columnar, and the strong 
repulsion phases. Furthermore, we checked the level- 
splitting conditions and evaluated the value of the cen- 
tral charge, which provided solid evidence to confirm the 
universality of the phase transition. By contrast, for the 
repulsive case, although we determined the dimer-liquid 
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APPENDIX A: A TRANSFER-MATRIX 
CALCULATION IN ONE-STRING SECTOR 



In this appendix, we shall explain how to construct the 
transfer matrix in the one-string sector and an analytical 



calculation of eigenvalues by the use of a Fourier transfor- 
mation. Since a row of the reference configuration given 
in Fig.[T]Jb) expresses a string vacuum state, we write it as 
|0) (see the top left in Fig. [8]). For convenience, we treat 
two sites in the x direction as one unit cell which includes 
four bonds (see squares by dotted blue lines). Then, one- 
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string states can be obtained via replacements of one of 
L c (= L/2) unit cells in |0) by four possible dimer con- 
figurations. From left to right of the first line in Fig. [8] 
(except for the vacuum), we call these as \Ai), \Bi), |C/), 
and \D{), respectively. Here, the center is supposed to 
be an Ith unit cell (I G [1,L C ]). Now, consider transfers 
of one-string states to those in the next-nearest-neighbor 
row in the y direction. Then, one can find 24 microscopic 



processes, which are listed in subsequent lines in Fig. [8] 
For instance, the second line shows seven microscopic 
processes of transfers from \Ai) in the first row to states 
in the third row, and thus exhibits an operation of the 
transfer matrix, i.e., T 2 _ L \A\). Consequently, we can 
obtain the following recursion relations for the transfers 
in the one-string sector 



J 



Tf_ Lc \At 



v 2 \A{) + h \B t ) + h \d) + |A) + \B l+1 ) + 



|CU) 

T'f_ Lc \Bi) =h \Q-i) + v 2 \Ai) + h 2 \Bi) + h \Ci) + h |A) + \B l+l ) + | D l+l ) , 
Ti_ Lc \Q) = + v 2 \Ai) + h \Bi) + h 2 \Q) + | A) + h \B l+1 ) + h \D l+1 ) , 

T?_ io | A) =v 2 |C,-i) + v 2 \Bi) + v 2 | A) , 



(Al) 
(A2) 
(A3) 
(A4) 



r 



where coefficients represent the Boltzmann weights of in- 
teractions in the first and the second rows. Next, by 
the use of the Fourier transformation, we can block- 
diagonalize the representation of Tf_ L ■ suppose that 



\X q ) = 4i=fl e ~ %ql \ Xl ) {X = A,B,C,D), (A5) 
v L < i=i 

then the q-block representation spanned by states 

\B q ), | Co), |Aj)} is given by a 4 x 4 complex nonsym- 



metric matrix: 



(v 2 



v 2 h 2 + z, 



v 

Vo 



g h + hz q 



h + hz q h + Zq 



l + z q 

h + Z q 

1 + hz 



(A6) 



■qi 



= e ±lq . Hence, the characteristic equation 



where 

to determine eigenvalues p is given by 



p 4 - 2(cosg + v 2 + h 2 )p 3 + [v 4 - 4h(l - h)v 2 + (1 - h 2 ) 2 ]p 2 + 2v 2 (l - h) 2 (l - v 2 - h 2 )p + w 4 (l - /i) 4 = 0. (A7) 



r 



Since it is invariant under a transformation q — > —q, a q 
dependence of the eigenvalue structure is even with re- 
spect to the point q = 0. Meanwhile, in general cases 
we use a software to evaluate q dependences of eigen- 
values; in some limiting cases, Eq. (|A7j) becomes simple 
and permits us to easily manipulate: for instance, for the 
noninteracting case h = v = 1, two of the four eigenval- 
ues are zero, and the rest is obtained from an equation 



p 2 — 2(cosq + 2)p +1 = 0. It then provides two real 
bands, as given in Fig. [7ta). In contrast, for the strong 
repulsion limit h — v = 0, two of the four eigenvalues are 
zero again, but others are complex values with a modu- 
lus of 1, i.e., e ±tq [see Fig. EJd)]. An implication of this 
eigenvalue structure, in particular a correspondence to 
the degeneracy of states in the IDM is discussed in the 
last part of Sec. IIIII 
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